Purpose
This document will analyze and display the stastistics of diversity at the locus and genotype level.
Required packages and data
library(PramCurry)
## Loading required package: poppr
## Loading required package: adegenet
## Loading required package: ade4
## ==========================
## adegenet 1.4-2 is loaded
## ==========================
##
## - to start, type '?adegenet'
## - to browse adegenet website, type 'adegenetWeb()'
## - to post questions/comments: adegenet-forum@lists.r-forge.r-project.org
##
##
## This is poppr version 1.1.2.99.70. To get started, type package?poppr
library(reshape2)
library(ggplot2)
library(poppr)
library(adegenet)
library(mmod)
## Loading required package: pegas
## Loading required package: ape
##
## Attaching package: 'pegas'
##
## The following object is masked from 'package:ape':
##
## mst
##
## The following object is masked from 'package:ade4':
##
## amova
library(dplyr)
##
## Attaching package: 'dplyr'
##
## The following objects are masked from 'package:stats':
##
## filter, lag
##
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(tidyr)
options(stringsAsFactors = FALSE)
data(ramdat)
data(pop_data)
data(myPal)
sessionInfo()
## R version 3.1.2 (2014-10-31)
## Platform: x86_64-apple-darwin13.4.0 (64-bit)
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] tidyr_0.1 dplyr_0.2 mmod_1.2.1
## [4] pegas_0.6 ape_3.1-4.5 ggplot2_1.0.0
## [7] reshape2_1.4 PramCurry_1.0 poppr_1.1.2.99-70
## [10] adegenet_1.4-2 ade4_1.6-2
##
## loaded via a namespace (and not attached):
## [1] assertthat_0.1 bitops_1.0-6 boot_1.3-11
## [4] caTools_1.17.1 colorspace_1.2-4 digest_0.6.4
## [7] evaluate_0.5.5 fastmatch_1.0-4 formatR_1.0
## [10] gdata_2.13.3 ggmap_2.3 grid_3.1.2
## [13] gtable_0.1.2 gtools_3.4.1 htmltools_0.2.6
## [16] httpuv_1.3.0 igraph_0.7.1 knitr_1.6
## [19] labtools_0.4.1 lattice_0.20-29 lubridate_1.3.3
## [22] mapproj_1.2-2 maps_2.3-9 MASS_7.3-34
## [25] Matrix_1.1-4 memoise_0.2.1 munsell_0.4.2
## [28] nlme_3.1-117 parallel_3.1.2 permute_0.8-3
## [31] phangorn_1.99-7 plotrix_3.5-7 plyr_1.8.1
## [34] png_0.1-7 proto_0.3-10 Rcpp_0.11.2
## [37] RgoogleMaps_1.2.0.6 rjson_0.2.14 RJSONIO_1.3-0
## [40] rmarkdown_0.3.10 scales_0.2.4 shiny_0.10.1
## [43] stringr_0.6.2 tools_3.1.2 vegan_2.0-10
## [46] xtable_1.7-4 yaml_2.1.13
myTheme <- theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5))
Functions
get_stats <- function(z){
mat <- matrix(nrow = nrow(z), ncol = 4,
dimnames = list(Pop = rownames(z),
Index = c("H", "G", "Hexp", "E.5")
)
)
N <- rowSums(z)
H <- vegan::diversity(z)
G <- vegan::diversity(z, "inv")
Simp <- vegan::diversity(z, "simp")
nei <- (N/(N-1)) * Simp
E.5 <- (G - 1)/(exp(H) -1)
mat[] <- c(H, G, nei, E.5)
return(mat)
}
boot_stats <- function(x, i){
res <- numeric(4)
names(res) <- c("H", "G", "Hexp", "E.5")
z <- table(x[i])
N <- sum(z)
H <- vegan::diversity(z)
G <- vegan::diversity(z, "inv")
Simp <- vegan::diversity(z, "simp")
nei <- (N/(N-1)) * Simp
E.5 <- (G - 1)/(exp(H) -1)
res[] <- c(H, G, nei, E.5)
return(res)
}
extract_samples <- function(x) rep(1:length(x), x)
do_boot <- function(tab, n, ...){
res <- apply(tab, 1, function(x) boot::boot(extract_samples(x), boot_stats, n, ...))
return(res)
}
get_ci <- function(x, lb, ub){
res <- apply(x$t, 2, quantile, c(lb, ub), na.rm = TRUE)
return(res)
}
get_all_ci <- function(res, ci = 95){
lower_bound <- (100 - ci)/200
upper_bound <- 1 - lower_bound
funval <- matrix(numeric(8), nrow = 2)
CI <- vapply(res, FUN = get_ci, FUN.VALUE = funval,
lower_bound, upper_bound)
dCI <- dimnames(CI)
dimnames(CI) <- list(CI = dCI[[1]],
Index = c("H", "G", "Hexp", "E.5"),
Pop = dCI[[3]])
return(CI)
}
boot_ci <- function(tab, n = 1000, ci = 95, total = TRUE, ...){
if (!is.matrix(tab) & is.genind(tab)){
tab <- mlg.table(tab, total = total, bar = FALSE)
}
res <- do_boot(tab, n, ...)
orig <- get_stats(tab)
orig <- melt(orig)
CI <- get_all_ci(res, ci = ci)
samp <- vapply(res, "[[", FUN.VALUE = res[[1]]$t, "t")
dimnames(samp) <- list(NULL,
Index = c("H", "G", "Hexp", "E.5"),
Pop = rownames(tab))
pl <- ggplot(melt(samp), aes(x = factor(Pop), y = value, group = Pop)) +
geom_boxplot() +
geom_point(aes(color = factor(Pop), x = factor(Pop), y = value),
size = 5, pch = 16, data = orig) +
xlab("Population") + labs(color = "Observed") +
facet_wrap(~Index, scales = "free_y") + myTheme
print(pl)
return(CI)
}
calc_loc_table <- function(dat, hier, combine = FALSE, lev = "allele", cc = TRUE){
if (cc){
dat <- clonecorrect(dat, hier = hier, combine = combine)
}
theTable <- locus_table(dat, information = FALSE, lev = lev)
theHier <- gethierarchy(dat, hier, combine = combine)
if (combine){
setpop(dat) <- hier
theHier <- levels(theHier[[length(theHier)]])
} else {
setpop(dat) <- as.formula(paste0("~", names(theHier)[1]))
theHier <- levels(theHier[[1]])
}
loc_tab <- lapply(theHier,
function(x){
locus_table(dat,
population = x,
information = FALSE,
lev = lev)
})
names(loc_tab) <- theHier
return(list(pops = loc_tab, total = theTable))
}
plot_loc_table <- function(tab){
# tab <- theList[["pops"]]
# tot <- melt(theList["total"])
tab <- melt(tab)
ggplot(tab, aes(y = value, x = L2, fill = L2, linetype = NA)) +
geom_bar(stat = "identity") +
geom_hline(aes(yintercept = value, linetype = L1),
data = tab[tab$L1 == "total", ]) +
facet_grid(summary ~ locus, scales = "free_y") +
scale_linetype_discrete(labels = "", name = "Pooled") +
myTheme
}
amova_pair <- function(pops, dat, hier, ...){
dat <- popsub(dat, pops)
if (any(table(pop(dat)) < 3)){
return(NULL)
}
poppr.amova(dat, hier, ...)
}
pairwise_amova <- function(x, hier, ...){
pops <- x@pop.names
pop_combs <- combn(pops, 2)
xlist <- apply(pop_combs, 2, amova_pair, x, hier, ...)
names(xlist) <- apply(pop_combs, 2, paste, collapse = " : ")
return(xlist)
}
pairwise_amova_test <- function(pairlist, nrepet = 99){
res <- lapply(names(pairlist), print_amova_pairs, pairlist, nrepet)
names(res) <- names(pairlist)
return(res)
}
print_amova_pairs <- function(pairname, pairlist, nrepet){
thePair <- pairlist[[pairname]]
if (!is.null(thePair)){
cat(pairname, "\n")
theTest <- randtest(thePair, nrepet = nrepet)
try(plot(theTest))
return(theTest)
} else {
return(NULL)
}
}
getpval <- function(x){
if (is.null(x$pvalue)){
return(as.numeric(rep(NA, 4)))
} else {
pvals <- x$pvalue
}
pvals[x$rep == 0] <- NA
return(pvals)
}
getphi <- function(x){
if (is.null(x)) return(numeric(4))
return(x$statphi$Phi)
}
refactor <- function(strings, factors){
u_factors <- unique(factors)
u_factors[match(strings, u_factors)]
}
make_amova_table <- function(am, amt){
tot <- nrow(am$results)
res <- data.frame(list(am$results[-tot, c("Df", "Sum Sq")],
Percent = am$componentsofcovariance[-tot, 2],
Pval = rev(amt$pvalue),
Phi = am$statphi$Phi[-tot]))
res <- as.matrix(res)
colnames(res) <- c("d.f.", "Sum of Squares", "Percent variation", "P",
"Phi statistic")
names(dimnames(res)) <- c("levels", "statistic")
return(res)
}
make_amova_printable <- function(amtab, amtabcc){
am_array <- array(dim = c(dim(amtab), 2),
dimnames = c(dimnames(amtab),
list(c("full", "clone-corrected"))))
am_array[, , 1] <- amtab
am_array[, , 2] <- amtabcc
tabfun <- function(x){
x <- paste0(paste0(signif(x, 3), collapse = " ("), ")")
return(x)
}
res <- apply(am_array, c(1, 2), tabfun)
return(res)
}
ci2string <- function(x, sig = 3, colString = "-"){
if (all(is.na(x))) return("(-)")
paste0("(", paste(signif(x, sig), collapse = colString),")")
}
add_ci <- function(dtable, ciarray, sig = 3, colString = "-"){
dtable <- data.frame(lapply(dtable, function(x, sig){
if (is.numeric(x)) x <- signif(x, sig)
return(x)
}, sig))
for (i in colnames(ciarray)){
to_add <- apply(ciarray[, i, ], 2, ci2string, sig, colString)
dtable[[i]] <- paste(signif(dtable[[i]], sig), to_add)
}
return(dtable)
}
Genotypic Diversity
(GDxYear <- poppr(ramdat,
sample = 999,
quiet = TRUE,
hist = FALSE))
## Pop N MLG eMLG SE H G Hexp E.5 Ia p.Ia rbarD
## 1 2001 39 10 3.65 1.161 1.187 1.90 0.486 0.395 -0.0864 0.869 -0.06765
## 2 2002 36 9 4.20 1.133 1.373 2.34 0.589 0.454 0.4555 0.001 0.19332
## 3 2003 102 20 4.75 1.301 1.811 2.92 0.665 0.376 0.4262 0.001 0.14326
## 4 2004 35 8 4.56 1.007 1.575 3.53 0.738 0.661 0.5282 0.001 0.18405
## 5 2005 2 2 2.00 0.000 0.693 2.00 1.000 1.000 NA NA NA
## 6 2006 1 1 1.00 0.000 0.000 1.00 NaN NaN NA NA NA
## 7 2010 1 1 1.00 0.000 0.000 1.00 NaN NaN NA NA NA
## 8 2011 68 6 2.15 0.827 0.563 1.32 0.245 0.420 0.1720 0.001 0.18094
## 9 2012 20 7 5.01 0.916 1.623 4.00 0.789 0.737 -0.1417 0.869 -0.17493
## 10 2013 179 47 7.74 1.204 3.148 13.77 0.933 0.573 0.0623 0.015 0.02925
## 11 2014 30 17 8.09 1.029 2.671 12.16 0.949 0.830 -0.0112 0.503 -0.00617
## 12 Total 513 70 6.95 1.332 2.983 8.64 0.886 0.408 0.2620 0.001 0.08159
## p.rD File
## 1 0.999 ramdat
## 2 0.001 ramdat
## 3 0.001 ramdat
## 4 0.001 ramdat
## 5 NA ramdat
## 6 NA ramdat
## 7 NA ramdat
## 8 0.001 ramdat
## 9 1.000 ramdat
## 10 0.001 ramdat
## 11 0.555 ramdat
## 12 0.001 ramdat
(GDxRegion <- poppr(setpop(ramdat, ~ZONE2),
sample = 999,
quiet = TRUE,
hist = FALSE))
## Pop N MLG eMLG SE H G Hexp E.5 Ia p.Ia
## 1 JHallCr 244 30 4.89 1.335 1.972 3.13 0.683 0.345 0.3322 0.001
## 2 NFChetHigh 114 35 6.81 1.342 2.735 7.07 0.866 0.421 0.1426 0.001
## 3 Coast 34 12 5.91 1.169 2.048 5.56 0.845 0.675 0.2514 0.011
## 4 HunterCr 66 4 1.88 0.690 0.423 1.24 0.198 0.460 -0.0440 0.770
## 5 Winchuck 35 9 4.29 1.072 1.474 2.88 0.672 0.559 -0.0196 0.593
## 6 ChetcoMain 16 7 5.00 0.926 1.450 2.84 0.692 0.565 0.3748 0.017
## 7 PistolRSF 4 2 2.00 0.000 0.562 1.60 0.500 0.795 0.0000 0.459
## 8 Total 513 70 6.95 1.332 2.983 8.64 0.886 0.408 0.2620 0.001
## rbarD p.rD File
## 1 0.1031 0.001 setpop(ramdat, ~ZONE2)
## 2 0.0660 0.001 setpop(ramdat, ~ZONE2)
## 3 0.2828 0.001 setpop(ramdat, ~ZONE2)
## 4 -0.0471 1.000 setpop(ramdat, ~ZONE2)
## 5 -0.0138 0.756 setpop(ramdat, ~ZONE2)
## 6 0.3957 0.001 setpop(ramdat, ~ZONE2)
## 7 NaN NA setpop(ramdat, ~ZONE2)
## 8 0.0816 0.001 setpop(ramdat, ~ZONE2)
(GDxYearccYR <- poppr(ramdat,
clonecorrect = TRUE,
hier = ~Pop/ZONE2,
sample = 999,
quiet = TRUE,
hist = FALSE))
## Pop N MLG eMLG SE H G Hexp E.5 Ia p.Ia rbarD
## 1 2001 10 10 10.00 0.000 2.303 10.0 1.000 1.000 -0.3622 0.961 -0.1891
## 2 2002 9 9 9.00 0.000 2.197 9.0 1.000 1.000 0.1957 0.207 0.0668
## 3 2003 21 20 9.79 0.410 2.979 19.2 0.995 0.974 0.2319 0.029 0.0615
## 4 2004 8 8 8.00 0.000 2.079 8.0 1.000 1.000 0.1611 0.262 0.0544
## 5 2005 2 2 2.00 0.000 0.693 2.0 1.000 1.000 NA NA NA
## 6 2006 1 1 1.00 0.000 0.000 1.0 NaN NaN NA NA NA
## 7 2010 1 1 1.00 0.000 0.000 1.0 NaN NaN NA NA NA
## 8 2011 6 6 6.00 0.000 1.792 6.0 1.000 1.000 -0.3023 0.783 -0.3226
## 9 2012 7 7 7.00 0.000 1.946 7.0 1.000 1.000 -0.3223 0.867 -0.3334
## 10 2013 68 47 9.42 0.710 3.724 35.6 0.986 0.855 -0.1275 0.998 -0.0508
## 11 2014 19 17 9.47 0.598 2.799 15.7 0.988 0.953 -0.1058 0.818 -0.0545
## 12 Total 152 70 9.19 0.829 3.939 38.4 0.980 0.742 0.0585 0.040 0.0155
## p.rD File
## 1 1.000 ramdat
## 2 0.145 ramdat
## 3 0.024 ramdat
## 4 0.212 ramdat
## 5 NA ramdat
## 6 NA ramdat
## 7 NA ramdat
## 8 1.000 ramdat
## 9 1.000 ramdat
## 10 1.000 ramdat
## 11 0.978 ramdat
## 12 0.035 ramdat
(GDxRegionccR <- poppr(ramdat,
clonecorrect = TRUE,
hier = ~ZONE2,
sample = 999,
quiet = TRUE,
hist = FALSE))
## Pop N MLG eMLG SE H G Hexp E.5 Ia p.Ia
## 1 JHallCr 30 30 10.0 NaN 3.401 30.0 1.000 1.000 0.0537 0.251
## 2 NFChetHigh 35 35 10.0 8.89e-07 3.555 35.0 1.000 1.000 -0.2156 1.000
## 3 Coast 12 12 10.0 NaN 2.485 12.0 1.000 1.000 0.0248 0.388
## 4 HunterCr 4 4 4.0 0.00e+00 1.386 4.0 1.000 1.000 -0.5714 0.855
## 5 Winchuck 9 9 9.0 0.00e+00 2.197 9.0 1.000 1.000 -0.3824 0.961
## 6 ChetcoMain 7 7 7.0 0.00e+00 1.946 7.0 1.000 1.000 -0.0945 0.581
## 7 PistolRSF 2 2 2.0 0.00e+00 0.693 2.0 1.000 1.000 NA NA
## 8 Total 99 70 9.6 6.06e-01 4.113 51.9 0.991 0.846 0.0409 0.146
## rbarD p.rD File
## 1 0.0136 0.250 ramdat
## 2 -0.0803 1.000 ramdat
## 3 0.0249 0.239 ramdat
## 4 -0.5774 0.999 ramdat
## 5 -0.1954 1.000 ramdat
## 6 -0.0988 0.924 ramdat
## 7 NA NA ramdat
## 8 0.0107 0.135 ramdat
(Totcc <- poppr(ramdat,
clonecorrect = TRUE,
hier = NA,
sublist = "Total",
sample = 999,
quiet = TRUE,
hist = FALSE))
## Pop N MLG eMLG SE H G Hexp E.5 Ia p.Ia rbarD p.rD
## 1 Total 70 70 70 0 4.25 70 1 1 -0.0216 0.659 -0.00551 0.663
## File
## 1 ramdat
Confidence intervals for genotypic diversity
# By Year:
set.seed(9001)
system.time(year_simp <- boot_ci(ramdat, n = 9999, total = TRUE,
parallel = "multicore", ncpus = 4L))
## Warning: Removed 19998 rows containing non-finite values (stat_boxplot).
## Warning: Removed 25017 rows containing non-finite values (stat_boxplot).
## Warning: Removed 2 rows containing missing values (geom_point).
## Warning: Removed 2 rows containing missing values (geom_point).

## user system elapsed
## 82.880 7.476 36.146
# By Population:
set.seed(9001)
system.time(zone_simp <- boot_ci(setpop(ramdat, ~ZONE2), n = 9999, total = TRUE,
parallel = "multicore", ncpus = 4L))
## Warning: Removed 3196 rows containing non-finite values (stat_boxplot).

## user system elapsed
## 61.808 6.224 26.320
write.table(add_ci(GDxYear, year_simp, colString = "-")[-nrow(GDxYear), -length(GDxYear)],
file = "GD.csv",
na = "-",
row.names = FALSE,
sep = ",", col.names = TRUE)
write.table(add_ci(GDxRegion, zone_simp, colString = "-")[-length(GDxYear)],
file = "GD.csv",
na = "-",
row.names = FALSE,
sep = ",", col.names = FALSE, append = TRUE)
Allelic and Genotypic diversity per locus
(loc_whole <- locus_table(ramdat)) # Whole data set
##
## allele = Number of observed alleles
## 1-D = Simpson index
## Hexp = Nei's 1978 expected heterozygosity
## ------------------------------------------
## summary
## locus allele 1-D Hexp Evenness
## PrMS6A1 2.00 0.50 1.00 1.00
## Pr9C3A1 2.00 0.50 1.00 1.00
## PrMS39A1 5.00 0.60 0.75 0.84
## PrMS45A1 3.00 0.50 0.75 0.99
## PrMS43A1 18.00 0.81 0.86 0.67
## mean 6.00 0.58 0.87 0.90
(loc_geno <- locus_table(clonecorrect(ramdat, hier = NA))) # By genotype
##
## allele = Number of observed alleles
## 1-D = Simpson index
## Hexp = Nei's 1978 expected heterozygosity
## ------------------------------------------
## summary
## locus allele 1-D Hexp Evenness
## PrMS6A1 2.00 0.50 0.99 0.99
## Pr9C3A1 2.00 0.49 0.99 0.99
## PrMS39A1 5.00 0.62 0.78 0.78
## PrMS45A1 3.00 0.51 0.76 0.95
## PrMS43A1 18.00 0.89 0.95 0.78
## mean 6.00 0.60 0.89 0.90
loc_year_allele <- calc_loc_table(ramdat, ~Pop/ZONE2)
(lya <- plot_loc_table(loc_year_allele) +
ggtitle("Allelic Diversity") + labs(list(x = "Year", fill = "Year")))
## Warning: Removed 1 rows containing missing values (position_stack).
## Warning: Removed 1 rows containing missing values (position_stack).
## Warning: Removed 1 rows containing missing values (position_stack).
## Warning: Removed 1 rows containing missing values (position_stack).
## Warning: Removed 1 rows containing missing values (position_stack).
## Warning: Removed 1 rows containing missing values (position_stack).
## Warning: Removed 1 rows containing missing values (position_stack).
## Warning: Removed 1 rows containing missing values (position_stack).
## Warning: Removed 1 rows containing missing values (position_stack).
## Warning: Removed 1 rows containing missing values (position_stack).
## Warning: Removed 1 rows containing missing values (position_stack).
## Warning: Removed 1 rows containing missing values (position_stack).
## Warning: Removed 1 rows containing missing values (position_stack).
## Warning: Removed 1 rows containing missing values (position_stack).
## Warning: Removed 1 rows containing missing values (position_stack).
## Warning: Removed 1 rows containing missing values (position_stack).
## Warning: Removed 1 rows containing missing values (position_stack).
## Warning: Removed 1 rows containing missing values (position_stack).
## Warning: Removed 1 rows containing missing values (position_stack).
## Warning: Removed 1 rows containing missing values (position_stack).
## Warning: Removed 1 rows containing missing values (position_stack).
## Warning: Removed 1 rows containing missing values (position_stack).
## Warning: Removed 1 rows containing missing values (position_stack).
## Warning: Removed 1 rows containing missing values (position_stack).

loc_year_geno <- calc_loc_table(ramdat, ~Pop/ZONE2, lev = "genotype")
(lyg <- plot_loc_table(loc_year_geno) +
ggtitle("Genotype Diversity") + labs(list(x = "Year", fill = "Year")))
## Warning: Removed 1 rows containing missing values (position_stack).
## Warning: Removed 1 rows containing missing values (position_stack).
## Warning: Removed 1 rows containing missing values (position_stack).
## Warning: Removed 1 rows containing missing values (position_stack).
## Warning: Removed 1 rows containing missing values (position_stack).
## Warning: Removed 1 rows containing missing values (position_stack).
## Warning: Removed 1 rows containing missing values (position_stack).
## Warning: Removed 1 rows containing missing values (position_stack).
## Warning: Removed 1 rows containing missing values (position_stack).
## Warning: Removed 1 rows containing missing values (position_stack).
## Warning: Removed 1 rows containing missing values (position_stack).
## Warning: Removed 1 rows containing missing values (position_stack).
## Warning: Removed 7 rows containing missing values (position_stack).
## Warning: Removed 8 rows containing missing values (position_stack).
## Warning: Removed 4 rows containing missing values (position_stack).
## Warning: Removed 8 rows containing missing values (position_stack).
## Warning: Removed 3 rows containing missing values (position_stack).
## Warning: Removed 3 rows containing missing values (position_stack).
## Warning: Removed 7 rows containing missing values (position_stack).
## Warning: Removed 8 rows containing missing values (position_stack).
## Warning: Removed 4 rows containing missing values (position_stack).
## Warning: Removed 8 rows containing missing values (position_stack).
## Warning: Removed 3 rows containing missing values (position_stack).
## Warning: Removed 3 rows containing missing values (position_stack).

loc_region_allele <- calc_loc_table(ramdat, ~ZONE2)
(lya <- plot_loc_table(loc_region_allele) +
ggtitle("Allelic Diversity") + labs(list(x = "Region", fill = "Region")))
## Warning: Removed 1 rows containing missing values (position_stack).
## Warning: Removed 1 rows containing missing values (position_stack).
## Warning: Removed 1 rows containing missing values (position_stack).
## Warning: Removed 1 rows containing missing values (position_stack).
## Warning: Removed 1 rows containing missing values (position_stack).
## Warning: Removed 1 rows containing missing values (position_stack).
## Warning: Removed 1 rows containing missing values (position_stack).
## Warning: Removed 1 rows containing missing values (position_stack).
## Warning: Removed 1 rows containing missing values (position_stack).
## Warning: Removed 1 rows containing missing values (position_stack).
## Warning: Removed 1 rows containing missing values (position_stack).
## Warning: Removed 1 rows containing missing values (position_stack).
## Warning: Removed 1 rows containing missing values (position_stack).
## Warning: Removed 1 rows containing missing values (position_stack).
## Warning: Removed 1 rows containing missing values (position_stack).
## Warning: Removed 1 rows containing missing values (position_stack).
## Warning: Removed 1 rows containing missing values (position_stack).
## Warning: Removed 1 rows containing missing values (position_stack).
## Warning: Removed 1 rows containing missing values (position_stack).
## Warning: Removed 1 rows containing missing values (position_stack).
## Warning: Removed 1 rows containing missing values (position_stack).
## Warning: Removed 1 rows containing missing values (position_stack).
## Warning: Removed 1 rows containing missing values (position_stack).
## Warning: Removed 1 rows containing missing values (position_stack).

loc_region_geno <- calc_loc_table(ramdat, ~ZONE2, lev = "genotype")
(lyg <- plot_loc_table(loc_region_geno) +
ggtitle("Genotype Diversity") + labs(list(x = "Region", fill = "Region")))
## Warning: Removed 1 rows containing missing values (position_stack).
## Warning: Removed 1 rows containing missing values (position_stack).
## Warning: Removed 1 rows containing missing values (position_stack).
## Warning: Removed 1 rows containing missing values (position_stack).
## Warning: Removed 1 rows containing missing values (position_stack).
## Warning: Removed 1 rows containing missing values (position_stack).
## Warning: Removed 1 rows containing missing values (position_stack).
## Warning: Removed 1 rows containing missing values (position_stack).
## Warning: Removed 1 rows containing missing values (position_stack).
## Warning: Removed 1 rows containing missing values (position_stack).
## Warning: Removed 1 rows containing missing values (position_stack).
## Warning: Removed 1 rows containing missing values (position_stack).
## Warning: Removed 6 rows containing missing values (position_stack).
## Warning: Removed 7 rows containing missing values (position_stack).
## Warning: Removed 2 rows containing missing values (position_stack).
## Warning: Removed 5 rows containing missing values (position_stack).
## Warning: Removed 1 rows containing missing values (position_stack).
## Warning: Removed 1 rows containing missing values (position_stack).
## Warning: Removed 6 rows containing missing values (position_stack).
## Warning: Removed 7 rows containing missing values (position_stack).
## Warning: Removed 2 rows containing missing values (position_stack).
## Warning: Removed 5 rows containing missing values (position_stack).
## Warning: Removed 1 rows containing missing values (position_stack).
## Warning: Removed 1 rows containing missing values (position_stack).

loc_region_allele_summary <- melt(loc_region_allele) %>%
filter(locus == "mean") %>%
mutate(region = ifelse(L1 == "pops", L2, L1)) %>%
select(locus, summary, value, region) %>%
dcast(region ~ summary) %>%
arrange(region)
loc_year_allele_summary <- melt(loc_year_allele) %>%
filter(locus == "mean") %>%
mutate(year = ifelse(L1 == "pops", L2, L1)) %>%
select(locus, summary, value, year) %>%
dcast(year ~ summary)
Tables
Mean Allelic Diversity by Region
| ChetcoMain |
3.0 |
0.5837 |
0.9655 |
0.9432 |
| Coast |
3.8 |
0.5944 |
0.9733 |
0.9597 |
| HunterCr |
2.6 |
0.5562 |
0.9615 |
0.9526 |
| JHallCr |
4.4 |
0.5720 |
0.9145 |
0.8908 |
| NFChetHigh |
4.8 |
0.5984 |
0.8961 |
0.8940 |
| PistolRSF |
2.4 |
0.5500 |
1.0000 |
1.0000 |
| Winchuck |
3.4 |
0.5630 |
0.9387 |
0.9165 |
| total |
6.0 |
0.6000 |
0.8909 |
0.9008 |
Mean Allelic Diversity by Year
| 2001 |
3.6 |
0.5640 |
0.9370 |
0.9099 |
| 2002 |
3.0 |
0.5642 |
0.9443 |
0.9306 |
| 2003 |
4.2 |
0.5646 |
0.9025 |
0.8931 |
| 2004 |
2.8 |
0.5031 |
0.8441 |
0.8548 |
| 2005 |
2.4 |
0.5000 |
0.9000 |
0.9180 |
| 2006 |
2.0 |
0.5000 |
1.0000 |
1.0000 |
| 2010 |
2.0 |
0.5000 |
1.0000 |
1.0000 |
| 2011 |
3.0 |
0.5833 |
0.9767 |
0.9638 |
| 2012 |
3.2 |
0.5755 |
0.9588 |
0.9498 |
| 2013 |
5.4 |
0.5979 |
0.8918 |
0.8987 |
| 2014 |
3.8 |
0.5909 |
0.9703 |
0.9347 |
| total |
6.0 |
0.5929 |
0.8823 |
0.8976 |
Total Allelic Diversity
| PrMS6A1 |
2.0000 |
0.4963 |
0.9927 |
0.9927 |
| Pr9C3A1 |
2.0000 |
0.4950 |
0.9900 |
0.9901 |
| PrMS39A1 |
5.0000 |
0.6229 |
0.7786 |
0.7817 |
| PrMS45A1 |
3.0000 |
0.5058 |
0.7587 |
0.9533 |
| PrMS43A1 |
18.0000 |
0.8932 |
0.9457 |
0.7837 |
| mean |
6.0000 |
0.6026 |
0.8931 |
0.9003 |
AMOVA
Including within individual variance
Year/Region
namehierarchy(ramdat) <- ~Year/ZONE2/ZONE1
(year_z2 <- poppr.amova(ramdat, ~Year/ZONE2))
##
## No missing values detected.
## $call
## ade4::amova(samples = xtab, distances = xdist, structures = xstruct)
##
## $results
## Df Sum Sq Mean Sq
## Between Year 10 79.79508 7.9795082
## Between ZONE2 Within Year 12 29.77109 2.4809238
## Between samples Within ZONE2 490 141.01083 0.2877772
## Within samples 513 1243.50000 2.4239766
## Total 1025 1494.07700 1.4576361
##
## $componentsofcovariance
## Sigma %
## Variations Between Year 0.02997373 2.035906
## Variations Between ZONE2 Within Year 0.08640422 5.868836
## Variations Between samples Within ZONE2 -1.06809970 -72.548560
## Variations Within samples 2.42397661 164.643817
## Total variations 1.47225486 100.000000
##
## $statphi
## Phi
## Phi-samples-total -0.64643817
## Phi-samples-ZONE2 -0.78775566
## Phi-ZONE2-Year 0.05990803
## Phi-Year-total 0.02035906
(year_z2cc <- poppr.amova(ramdat, ~Year/ZONE2, clonecorrect = TRUE))
##
## No missing values detected.
## $call
## ade4::amova(samples = xtab, distances = xdist, structures = xstruct)
##
## $results
## Df Sum Sq Mean Sq
## Between Year 10 10.523948 1.0523948
## Between ZONE2 Within Year 12 9.459087 0.7882573
## Between samples Within ZONE2 129 71.609070 0.5551091
## Within samples 152 359.000000 2.3618421
## Total 303 450.592105 1.4871027
##
## $componentsofcovariance
## Sigma %
## Variations Between Year 0.009753153 0.6553138
## Variations Between ZONE2 Within Year 0.020089101 1.3497856
## Variations Between samples Within ZONE2 -0.903366518 -60.6971504
## Variations Within samples 2.361842105 158.6920509
## Total variations 1.488317840 100.0000000
##
## $statphi
## Phi
## Phi-samples-total -0.586920509
## Phi-samples-ZONE2 -0.619390908
## Phi-ZONE2-Year 0.013586893
## Phi-Year-total 0.006553138
set.seed(9001)
system.time(year_z2.t <- randtest(year_z2, nrepet = 999))
## user system elapsed
## 153.881 1.489 155.397
plot(year_z2.t)

set.seed(9001)
system.time(year_z2cc.t <- randtest(year_z2cc, nrepet = 999))
## user system elapsed
## 17.137 0.212 17.353
plot(year_z2cc.t)

date()
## [1] "Fri Oct 24 16:21:38 2014"
Region/Year
namehierarchy(ramdat) <- ~Year/ZONE2/ZONE1
(z2_year <- poppr.amova(ramdat, ~ZONE2/Year))
##
## No missing values detected.
## $call
## ade4::amova(samples = xtab, distances = xdist, structures = xstruct)
##
## $results
## Df Sum Sq Mean Sq
## Between ZONE2 6 98.53739 16.4228985
## Between Year Within ZONE2 16 11.02878 0.6892985
## Between samples Within Year 490 141.01083 0.2877772
## Within samples 513 1243.50000 2.4239766
## Total 1025 1494.07700 1.4576361
##
## $componentsofcovariance
## Sigma %
## Variations Between ZONE2 0.1283283 8.5740970
## Variations Between Year Within ZONE2 0.0124923 0.8346573
## Variations Between samples Within Year -1.0680997 -71.3637662
## Variations Within samples 2.4239766 161.9550118
## Total variations 1.4966975 100.0000000
##
## $statphi
## Phi
## Phi-samples-total -0.619550118
## Phi-samples-Year -0.787755656
## Phi-Year-ZONE2 0.009129331
## Phi-ZONE2-total 0.085740970
(z2_yearcc <- poppr.amova(ramdat, ~ZONE2/Year, clonecorrect = TRUE))
##
## No missing values detected.
## $call
## ade4::amova(samples = xtab, distances = xdist, structures = xstruct)
##
## $results
## Df Sum Sq Mean Sq
## Between ZONE2 6 11.582493 1.9304155
## Between Year Within ZONE2 16 8.400542 0.5250339
## Between samples Within Year 129 71.609070 0.5551091
## Within samples 152 359.000000 2.3618421
## Total 303 450.592105 1.4871027
##
## $componentsofcovariance
## Sigma %
## Variations Between ZONE2 0.039160694 2.6194612
## Variations Between Year Within ZONE2 -0.002645858 -0.1769816
## Variations Between samples Within Year -0.903366518 -60.4262411
## Variations Within samples 2.361842105 157.9837615
## Total variations 1.494990423 100.0000000
##
## $statphi
## Phi
## Phi-samples-total -0.579837615
## Phi-samples-Year -0.619390908
## Phi-Year-ZONE2 -0.001817423
## Phi-ZONE2-total 0.026194612
set.seed(9001)
system.time(z2_year.t <- randtest(z2_year, nrepet = 999))
## user system elapsed
## 154.778 0.348 155.143
plot(z2_year.t)

set.seed(9001)
system.time(z2_yearcc.t <- randtest(z2_yearcc, nrepet = 999))
## user system elapsed
## 17.025 0.012 17.034
plot(z2_yearcc.t)

date()
## [1] "Fri Oct 24 16:42:51 2014"
Excluding within individual variance
Year/Region
namehierarchy(ramdat) <- ~Year/ZONE2/ZONE1
(WF_year_z2 <- poppr.amova(ramdat, ~Year/ZONE2, within = FALSE))
##
## No missing values detected.
## $call
## ade4::amova(samples = xtab, distances = xdist, structures = xstruct)
##
## $results
## Df Sum Sq Mean Sq
## Between Year 10 159.7 15.9691
## Between samples Within Year 12 59.5 4.9581
## Within samples 490 280.8 0.5731
## Total 512 500.0 0.9766
##
## $componentsofcovariance
## Sigma %
## Variations Between Year 0.1203 11.58
## Variations Between samples Within Year 0.3455 33.26
## Variations Within samples 0.5731 55.16
## Total variations 1.0389 100.00
##
## $statphi
## Phi
## Phi-samples-total 0.4484
## Phi-samples-Year 0.3761
## Phi-Year-total 0.1158
(WF_year_z2cc <- poppr.amova(ramdat, ~Year/ZONE2, clonecorrect = TRUE,
within = FALSE))
##
## No missing values detected.
## $call
## ade4::amova(samples = xtab, distances = xdist, structures = xstruct)
##
## $results
## Df Sum Sq Mean Sq
## Between Year 10 21.01 2.101
## Between samples Within Year 12 19.14 1.595
## Within samples 129 141.29 1.095
## Total 151 181.44 1.202
##
## $componentsofcovariance
## Sigma %
## Variations Between Year 0.03654 3.001
## Variations Between samples Within Year 0.08615 7.073
## Variations Within samples 1.09524 89.926
## Total variations 1.21794 100.000
##
## $statphi
## Phi
## Phi-samples-total 0.10074
## Phi-samples-Year 0.07292
## Phi-Year-total 0.03001
set.seed(9001)
system.time(WF_year_z2.t <- randtest(WF_year_z2, nrepet = 9999))
## user system elapsed
## 35.34 0.27 35.98
plot(WF_year_z2.t)

set.seed(9001)
system.time(WF_year_z2cc.t <- randtest(WF_year_z2cc, nrepet = 9999))
## user system elapsed
## 18.973 0.097 19.129
plot(WF_year_z2cc.t)

date()
## [1] "Tue Nov 4 09:32:42 2014"
Region/Year
namehierarchy(ramdat) <- ~Year/ZONE2/ZONE1
(WF_z2_year <- poppr.amova(ramdat, ~ZONE2/Year, within = FALSE))
##
## No missing values detected.
## $call
## ade4::amova(samples = xtab, distances = xdist, structures = xstruct)
##
## $results
## Df Sum Sq Mean Sq
## Between ZONE2 6 196.69 32.7822
## Between samples Within ZONE2 16 22.49 1.4059
## Within samples 490 280.82 0.5731
## Total 512 500.00 0.9766
##
## $componentsofcovariance
## Sigma %
## Variations Between ZONE2 0.51127 44.999
## Variations Between samples Within ZONE2 0.05182 4.561
## Variations Within samples 0.57309 50.440
## Total variations 1.13618 100.000
##
## $statphi
## Phi
## Phi-samples-total 0.49560
## Phi-samples-ZONE2 0.08293
## Phi-ZONE2-total 0.44999
(WF_z2_yearcc <- poppr.amova(ramdat, ~ZONE2/Year, clonecorrect = TRUE,
within = FALSE))
##
## No missing values detected.
## $call
## ade4::amova(samples = xtab, distances = xdist, structures = xstruct)
##
## $results
## Df Sum Sq Mean Sq
## Between ZONE2 6 22.95 3.825
## Between samples Within ZONE2 16 17.20 1.075
## Within samples 129 141.29 1.095
## Total 151 181.44 1.202
##
## $componentsofcovariance
## Sigma %
## Variations Between ZONE2 0.152442 12.2526
## Variations Between samples Within ZONE2 -0.003527 -0.2835
## Variations Within samples 1.095242 88.0308
## Total variations 1.244157 100.0000
##
## $statphi
## Phi
## Phi-samples-total 0.119692
## Phi-samples-ZONE2 -0.003231
## Phi-ZONE2-total 0.122526
set.seed(9001)
system.time(WF_z2_year.t <- randtest(WF_z2_year, nrepet = 9999))
## user system elapsed
## 35.077 0.341 35.820
plot(WF_z2_year.t)

set.seed(9001)
system.time(WF_z2_yearcc.t <- randtest(WF_z2_yearcc, nrepet = 9999))
## user system elapsed
## 17.879 0.099 18.057
plot(WF_z2_yearcc.t)

date()
## [1] "Tue Nov 4 09:33:37 2014"
Writing tables
z2year_table <- make_amova_table(WF_z2_year, WF_z2_year.t)
z2yearcc_table <- make_amova_table(WF_z2_yearcc, WF_z2_yearcc.t)
(z2year_full_table <- make_amova_printable(z2year_table, z2yearcc_table))
## statistic
## levels d.f. Sum of Squares
## Between ZONE2 "6 (6)" "197 (23)"
## Between samples Within ZONE2 "16 (16)" "22.5 (17.2)"
## Within samples "490 (129)" "281 (141)"
## statistic
## levels Percent variation P
## Between ZONE2 "45 (12.3)" "1e-04 (1e-04)"
## Between samples Within ZONE2 "4.56 (-0.283)" "1e-04 (0.446)"
## Within samples "50.4 (88)" "1e-04 (1e-04)"
## statistic
## levels Phi statistic
## Between ZONE2 "0.496 (0.12)"
## Between samples Within ZONE2 "0.0829 (-0.00323)"
## Within samples "0.45 (0.123)"
write.table(z2year_full_table, file = "zone_by_year.csv", row.names = TRUE,
col.names = NA, sep = ",")
yearz2_table <- make_amova_table(WF_year_z2, WF_year_z2.t)
yearz2cc_table <- make_amova_table(WF_year_z2cc, WF_year_z2cc.t)
(yearz2_full_table <- make_amova_printable(yearz2_table, yearz2cc_table))
## statistic
## levels d.f. Sum of Squares Percent variation
## Between Year "10 (10)" "160 (21)" "11.6 (3)"
## Between samples Within Year "12 (12)" "59.5 (19.1)" "33.3 (7.07)"
## Within samples "490 (129)" "281 (141)" "55.2 (89.9)"
## statistic
## levels P Phi statistic
## Between Year "0.366 (0.175)" "0.448 (0.101)"
## Between samples Within Year "1e-04 (2e-04)" "0.376 (0.0729)"
## Within samples "1e-04 (1e-04)" "0.116 (0.03)"
write.table(yearz2_full_table, file = "year_by_zone.csv", row.names = TRUE,
col.names = NA, sep = ",")
Without Cape Sebastian or Pistol River South Fork
noseb <- popsub(setpop(ramdat, ~ZONE2), blacklist = "HunterCr")
system.time(noseb_year_z2 <- poppr.amova(noseb, ~Year/ZONE2))
##
## No missing values detected.
## user system elapsed
## 0.776 0.000 0.778
system.time(noseb_year_z2cc <- poppr.amova(noseb, ~Year/ZONE2, clonecorrect = TRUE))
##
## No missing values detected.
## user system elapsed
## 0.892 0.000 0.892
set.seed(9001)
system.time(noseb_year_z2.t <- randtest(noseb_year_z2, nrepet = 999))
## user system elapsed
## 114.244 0.664 114.935
plot(noseb_year_z2.t)

set.seed(9001)
system.time(noseb_year_z2cc.t <- randtest(noseb_year_z2cc, nrepet = 999))
## user system elapsed
## 15.653 0.000 15.658
plot(noseb_year_z2cc.t)

date()
## [1] "Fri Oct 24 16:23:52 2014"
noseb <- popsub(setpop(ramdat, ~ZONE2), blacklist = "HunterCr")
system.time(noseb_z2_year <- poppr.amova(noseb, ~ZONE2/Year))
##
## No missing values detected.
## user system elapsed
## 0.768 0.052 0.821
system.time(noseb_z2_yearcc <- poppr.amova(noseb, ~ZONE2/Year, clonecorrect = TRUE))
##
## No missing values detected.
## user system elapsed
## 0.900 0.016 0.915
set.seed(9001)
system.time(noseb_z2_year.t <- randtest(noseb_z2_year, nrepet = 999))
## user system elapsed
## 113.627 1.380 115.031
plot(noseb_z2_year.t)

set.seed(9001)
system.time(noseb_z2_yearcc.t <- randtest(noseb_z2_yearcc, nrepet = 999))
## user system elapsed
## 15.425 0.200 15.629
plot(noseb_z2_yearcc.t)

date()
## [1] "Fri Oct 24 16:45:05 2014"
nosebpr <- selPopSize(noseb, n = 10)
system.time(nosebpr_year_z2 <- poppr.amova(nosebpr, ~Year/ZONE2))
##
## No missing values detected.
## user system elapsed
## 0.764 0.004 0.767
system.time(nosebpr_year_z2cc <- poppr.amova(nosebpr, ~Year/ZONE2, clonecorrect = TRUE))
##
## No missing values detected.
## user system elapsed
## 0.852 0.024 0.875
set.seed(9001)
system.time(nosebpr_year_z2.t <- randtest(nosebpr_year_z2, nrepet = 999))
## user system elapsed
## 112.171 1.472 113.661
plot(nosebpr_year_z2.t)

set.seed(9001)
system.time(nosebpr_year_z2cc.t <- randtest(nosebpr_year_z2cc, nrepet = 999))
## user system elapsed
## 15.269 0.216 15.489
plot(nosebpr_year_z2cc.t)

nosebpr <- selPopSize(noseb, n = 10)
system.time(nosebpr_z2_year <- poppr.amova(nosebpr, ~ZONE2/Year))
##
## No missing values detected.
## user system elapsed
## 0.780 0.012 0.793
system.time(nosebpr_z2_yearcc <- poppr.amova(nosebpr, ~ZONE2/Year, clonecorrect = TRUE))
##
## No missing values detected.
## user system elapsed
## 0.860 0.008 0.867
set.seed(9001)
system.time(nosebpr_z2_year.t <- randtest(nosebpr_z2_year, nrepet = 999))
## user system elapsed
## 110.799 1.376 112.190
plot(nosebpr_z2_year.t)

set.seed(9001)
system.time(nosebpr_z2_yearcc.t <- randtest(nosebpr_z2_yearcc, nrepet = 999))
## user system elapsed
## 15.165 0.160 15.326
plot(nosebpr_z2_yearcc.t)

Excluding individual variance and Cape Sebastian/PRSF
Year/Region
nosebpr <- resetMLG(nosebpr)
(nosebpr_WF_year_z2 <- poppr.amova(nosebpr, ~Year/ZONE2, within = FALSE))
##
## No missing values detected.
## $call
## ade4::amova(samples = xtab, distances = xdist, structures = xstruct)
##
## $results
## Df Sum Sq Mean Sq
## Between Year 10 48.86 4.8862
## Between samples Within Year 10 54.55 5.4547
## Within samples 422 272.72 0.6463
## Total 442 376.13 0.8510
##
## $componentsofcovariance
## Sigma %
## Variations Between Year -0.1154 -13.31
## Variations Between samples Within Year 0.3359 38.75
## Variations Within samples 0.6463 74.56
## Total variations 0.8668 100.00
##
## $statphi
## Phi
## Phi-samples-total 0.2544
## Phi-samples-Year 0.3420
## Phi-Year-total -0.1331
(nosebpr_WF_year_z2cc <- poppr.amova(nosebpr, ~Year/ZONE2, clonecorrect = TRUE,
within = FALSE))
##
## No missing values detected.
## $call
## ade4::amova(samples = xtab, distances = xdist, structures = xstruct)
##
## $results
## Df Sum Sq Mean Sq
## Between Year 10 18.95 1.895
## Between samples Within Year 10 16.47 1.647
## Within samples 125 138.04 1.104
## Total 145 173.45 1.196
##
## $componentsofcovariance
## Sigma %
## Variations Between Year 0.02078 1.718
## Variations Between samples Within Year 0.08416 6.960
## Variations Within samples 1.10429 91.322
## Total variations 1.20923 100.000
##
## $statphi
## Phi
## Phi-samples-total 0.08678
## Phi-samples-Year 0.07082
## Phi-Year-total 0.01718
set.seed(9001)
system.time(nosebpr_WF_year_z2.t <- randtest(nosebpr_WF_year_z2, nrepet = 9999))
## user system elapsed
## 27.845 0.185 28.338
plot(nosebpr_WF_year_z2.t)

set.seed(9001)
system.time(nosebpr_WF_year_z2cc.t <- randtest(nosebpr_WF_year_z2cc, nrepet = 9999))
## user system elapsed
## 16.948 0.205 19.099
plot(nosebpr_WF_year_z2cc.t)

date()
## [1] "Fri Nov 14 10:34:00 2014"
Region/Year
namehierarchy(nosebpr) <- ~Year/ZONE2/ZONE1
(nosebpr_WF_z2_year <- poppr.amova(nosebpr, ~ZONE2/Year, within = FALSE))
##
## No missing values detected.
## $call
## ade4::amova(samples = xtab, distances = xdist, structures = xstruct)
##
## $results
## Df Sum Sq Mean Sq
## Between ZONE2 4 80.91 20.2284
## Between samples Within ZONE2 16 22.49 1.4059
## Within samples 422 272.72 0.6463
## Total 442 376.13 0.8510
##
## $componentsofcovariance
## Sigma %
## Variations Between ZONE2 0.26437 27.599
## Variations Between samples Within ZONE2 0.04727 4.935
## Variations Within samples 0.64627 67.467
## Total variations 0.95791 100.000
##
## $statphi
## Phi
## Phi-samples-total 0.32533
## Phi-samples-ZONE2 0.06816
## Phi-ZONE2-total 0.27599
(nosebpr_WF_z2_yearcc <- poppr.amova(nosebpr, ~ZONE2/Year, clonecorrect = TRUE,
within = FALSE))
##
## No missing values detected.
## $call
## ade4::amova(samples = xtab, distances = xdist, structures = xstruct)
##
## $results
## Df Sum Sq Mean Sq
## Between ZONE2 4 18.21 4.553
## Between samples Within ZONE2 16 17.20 1.075
## Within samples 125 138.04 1.104
## Total 145 173.45 1.196
##
## $componentsofcovariance
## Sigma %
## Variations Between ZONE2 0.138799 11.2119
## Variations Between samples Within ZONE2 -0.005119 -0.4135
## Variations Within samples 1.104290 89.2016
## Total variations 1.237970 100.0000
##
## $statphi
## Phi
## Phi-samples-total 0.107984
## Phi-samples-ZONE2 -0.004657
## Phi-ZONE2-total 0.112119
set.seed(9001)
system.time(nosebpr_WF_z2_year.t <- randtest(nosebpr_WF_z2_year, nrepet = 9999))
## user system elapsed
## 26.774 0.203 27.382
plot(nosebpr_WF_z2_year.t)

set.seed(9001)
system.time(nosebpr_WF_z2_yearcc.t <- randtest(nosebpr_WF_z2_yearcc, nrepet = 9999))
## user system elapsed
## 14.217 0.091 14.394
plot(nosebpr_WF_z2_yearcc.t)

date()
## [1] "Fri Nov 14 10:34:43 2014"
nosebpr_z2year_table <- make_amova_table(nosebpr_WF_z2_year, nosebpr_WF_z2_year.t)
nosebpr_z2yearcc_table <- make_amova_table(nosebpr_WF_z2_yearcc, nosebpr_WF_z2_yearcc.t)
(nosebpr_z2year_full_table <- make_amova_printable(nosebpr_z2year_table, nosebpr_z2yearcc_table))
## statistic
## levels d.f. Sum of Squares
## Between ZONE2 "4 (4)" "80.9 (18.2)"
## Between samples Within ZONE2 "16 (16)" "22.5 (17.2)"
## Within samples "422 (125)" "273 (138)"
## statistic
## levels Percent variation P
## Between ZONE2 "27.6 (11.2)" "1e-04 (1e-04)"
## Between samples Within ZONE2 "4.93 (-0.413)" "1e-04 (0.445)"
## Within samples "67.5 (89.2)" "1e-04 (1e-04)"
## statistic
## levels Phi statistic
## Between ZONE2 "0.325 (0.108)"
## Between samples Within ZONE2 "0.0682 (-0.00466)"
## Within samples "0.276 (0.112)"
write.table(nosebpr_z2year_full_table, file = "nosebpr_zone_by_year.csv", row.names = TRUE,
col.names = NA, sep = ",")
nosebpr_yearz2_table <- make_amova_table(nosebpr_WF_year_z2, nosebpr_WF_year_z2.t)
nosebpr_yearz2cc_table <- make_amova_table(nosebpr_WF_year_z2cc, nosebpr_WF_year_z2cc.t)
(nosebpr_yearz2_full_table <- make_amova_printable(nosebpr_yearz2_table, nosebpr_yearz2cc_table))
## statistic
## levels d.f. Sum of Squares Percent variation
## Between Year "10 (10)" "48.9 (18.9)" "-13.3 (1.72)"
## Between samples Within Year "10 (10)" "54.5 (16.5)" "38.8 (6.96)"
## Within samples "422 (125)" "273 (138)" "74.6 (91.3)"
## statistic
## levels P Phi statistic
## Between Year "0.796 (0.293)" "0.254 (0.0868)"
## Between samples Within Year "1e-04 (4e-04)" "0.342 (0.0708)"
## Within samples "1e-04 (1e-04)" "-0.133 (0.0172)"
write.table(nosebpr_yearz2_full_table, file = "nosebpr_year_by_zone.csv", row.names = TRUE,
col.names = NA, sep = ",")
Without the coast, cape sebastian, or PRSF
Year/Region
nosebprc <- popsub(setpop(nosebpr, ~ZONE2), blacklist = "Coast")
nosebprc <- resetMLG(nosebprc)
(nosebprc_WF_year_z2 <- poppr.amova(nosebprc, ~Year/ZONE2, within = FALSE))
##
## No missing values detected.
## $call
## ade4::amova(samples = xtab, distances = xdist, structures = xstruct)
##
## $results
## Df Sum Sq Mean Sq
## Between Year 7 41.20 5.8852
## Between samples Within Year 7 37.96 5.4226
## Within samples 394 253.02 0.6422
## Total 408 332.17 0.8141
##
## $componentsofcovariance
## Sigma %
## Variations Between Year -0.1232 -14.83
## Variations Between samples Within Year 0.3116 37.52
## Variations Within samples 0.6422 77.31
## Total variations 0.8306 100.00
##
## $statphi
## Phi
## Phi-samples-total 0.2269
## Phi-samples-Year 0.3267
## Phi-Year-total -0.1483
(nosebprc_WF_year_z2cc <- poppr.amova(nosebprc, ~Year/ZONE2, clonecorrect = TRUE,
within = FALSE))
##
## No missing values detected.
## $call
## ade4::amova(samples = xtab, distances = xdist, structures = xstruct)
##
## $results
## Df Sum Sq Mean Sq
## Between Year 7 14.30 2.043
## Between samples Within Year 7 10.03 1.433
## Within samples 112 125.39 1.120
## Total 126 149.72 1.188
##
## $componentsofcovariance
## Sigma %
## Variations Between Year 0.03813 3.168
## Variations Between samples Within Year 0.04603 3.824
## Variations Within samples 1.11958 93.008
## Total variations 1.20375 100.000
##
## $statphi
## Phi
## Phi-samples-total 0.06992
## Phi-samples-Year 0.03949
## Phi-Year-total 0.03168
set.seed(9001)
system.time(nosebprc_WF_year_z2.t <- randtest(nosebprc_WF_year_z2, nrepet = 9999))
## user system elapsed
## 18.395 0.074 18.592
plot(nosebprc_WF_year_z2.t)

set.seed(9001)
system.time(nosebprc_WF_year_z2cc.t <- randtest(nosebprc_WF_year_z2cc, nrepet = 9999))
## user system elapsed
## 9.906 0.040 9.966
plot(nosebprc_WF_year_z2cc.t)

date()
## [1] "Sun Nov 16 09:51:22 2014"
Region/Year
namehierarchy(nosebprc) <- ~Year/ZONE2/ZONE1
(nosebprc_WF_z2_year <- poppr.amova(nosebprc, ~ZONE2/Year, within = FALSE))
##
## No missing values detected.
## $call
## ade4::amova(samples = xtab, distances = xdist, structures = xstruct)
##
## $results
## Df Sum Sq Mean Sq
## Between ZONE2 3 62.36 20.7878
## Between samples Within ZONE2 11 16.79 1.5265
## Within samples 394 253.02 0.6422
## Total 408 332.17 0.8141
##
## $componentsofcovariance
## Sigma %
## Variations Between ZONE2 0.24461 26.354
## Variations Between samples Within ZONE2 0.04136 4.456
## Variations Within samples 0.64217 69.189
## Total variations 0.92814 100.000
##
## $statphi
## Phi
## Phi-samples-total 0.30811
## Phi-samples-ZONE2 0.06051
## Phi-ZONE2-total 0.26354
(nosebprc_WF_z2_yearcc <- poppr.amova(nosebprc, ~ZONE2/Year, clonecorrect = TRUE,
within = FALSE))
##
## No missing values detected.
## $call
## ade4::amova(samples = xtab, distances = xdist, structures = xstruct)
##
## $results
## Df Sum Sq Mean Sq
## Between ZONE2 3 11.59 3.863
## Between samples Within ZONE2 11 12.74 1.158
## Within samples 112 125.39 1.120
## Total 126 149.72 1.188
##
## $componentsofcovariance
## Sigma %
## Variations Between ZONE2 0.103175 8.400
## Variations Between samples Within ZONE2 0.005491 0.447
## Variations Within samples 1.119584 91.153
## Total variations 1.228250 100.000
##
## $statphi
## Phi
## Phi-samples-total 0.08847
## Phi-samples-ZONE2 0.00488
## Phi-ZONE2-total 0.08400
set.seed(9001)
system.time(nosebprc_WF_z2_year.t <- randtest(nosebprc_WF_z2_year, nrepet = 9999))
## user system elapsed
## 17.590 0.057 17.679
plot(nosebprc_WF_z2_year.t)

set.seed(9001)
system.time(nosebprc_WF_z2_yearcc.t <- randtest(nosebprc_WF_z2_yearcc, nrepet = 9999))
## user system elapsed
## 9.205 0.036 9.258
plot(nosebprc_WF_z2_yearcc.t)

date()
## [1] "Sun Nov 16 09:51:53 2014"
nosebprc_z2year_table <- make_amova_table(nosebprc_WF_z2_year, nosebprc_WF_z2_year.t)
nosebprc_z2yearcc_table <- make_amova_table(nosebprc_WF_z2_yearcc, nosebprc_WF_z2_yearcc.t)
(nosebprc_z2year_full_table <- make_amova_printable(nosebprc_z2year_table, nosebprc_z2yearcc_table))
## statistic
## levels d.f. Sum of Squares
## Between ZONE2 "3 (3)" "62.4 (11.6)"
## Between samples Within ZONE2 "11 (11)" "16.8 (12.7)"
## Within samples "394 (112)" "253 (125)"
## statistic
## levels Percent variation P
## Between ZONE2 "26.4 (8.4)" "2e-04 (0.0012)"
## Between samples Within ZONE2 "4.46 (0.447)" "1e-04 (0.362)"
## Within samples "69.2 (91.2)" "1e-04 (1e-04)"
## statistic
## levels Phi statistic
## Between ZONE2 "0.308 (0.0885)"
## Between samples Within ZONE2 "0.0605 (0.00488)"
## Within samples "0.264 (0.084)"
write.table(nosebprc_z2year_full_table, file = "nosebprc_zone_by_year.csv", row.names = TRUE,
col.names = NA, sep = ",")
nosebprc_yearz2_table <- make_amova_table(nosebprc_WF_year_z2, nosebprc_WF_year_z2.t)
nosebprc_yearz2cc_table <- make_amova_table(nosebprc_WF_year_z2cc, nosebprc_WF_year_z2cc.t)
(nosebprc_yearz2_full_table <- make_amova_printable(nosebprc_yearz2_table, nosebprc_yearz2cc_table))
## statistic
## levels d.f. Sum of Squares Percent variation
## Between Year "7 (7)" "41.2 (14.3)" "-14.8 (3.17)"
## Between samples Within Year "7 (7)" "38 (10)" "37.5 (3.82)"
## Within samples "394 (112)" "253 (125)" "77.3 (93)"
## statistic
## levels P Phi statistic
## Between Year "0.843 (0.244)" "0.227 (0.0699)"
## Between samples Within Year "1e-04 (0.0265)" "0.327 (0.0395)"
## Within samples "1e-04 (3e-04)" "-0.148 (0.0317)"
write.table(nosebprc_yearz2_full_table, file = "nosebprc_year_by_zone.csv", row.names = TRUE,
col.names = NA, sep = ",")
Pairwise AMOVA Between Regions
While presenting the AMOVA results, Ebba brought up a very good point that the more recently introduced regions might be biasing the AMOVA results. To assess which regions might have differentiation in a pairwise fashion, a pairwise AMOVA is being run on regions comparing both regions separately per year.
Running the AMOVA by zone
zonePairs <- pairwise_amova(setpop(ramdat, ~ZONE2), ~ZONE2/Year, quiet = TRUE)
zonePairs.cc <- pairwise_amova(setpop(ramdat, ~ZONE2), ~ZONE2/Year,
quiet = TRUE, clonecorrect = TRUE)
set.seed(9001) #
system.time(zonePairs.t <- pairwise_amova_test(zonePairs, nrepet = 999))
## JHallCr : NFChetHigh

## JHallCr : Coast

## JHallCr : HunterCr

## JHallCr : Winchuck

## JHallCr : ChetcoMain

## JHallCr : PistolRSF

## NFChetHigh : Coast

## NFChetHigh : HunterCr

## NFChetHigh : Winchuck

## NFChetHigh : ChetcoMain

## NFChetHigh : PistolRSF

## Coast : HunterCr

## Coast : Winchuck

## Coast : ChetcoMain

## Coast : PistolRSF

## HunterCr : Winchuck

## HunterCr : ChetcoMain

## HunterCr : PistolRSF

## Winchuck : ChetcoMain

## Winchuck : PistolRSF

## ChetcoMain : PistolRSF

## user system elapsed
## 297.830 0.140 298.174
set.seed(9001)
system.time(zonePairs.cc.t <- pairwise_amova_test(zonePairs.cc, nrepet = 999))
## JHallCr : NFChetHigh

## JHallCr : Coast

## JHallCr : HunterCr

## JHallCr : Winchuck

## JHallCr : ChetcoMain

## JHallCr : PistolRSF

## NFChetHigh : Coast

## NFChetHigh : HunterCr

## NFChetHigh : Winchuck

## NFChetHigh : ChetcoMain

## NFChetHigh : PistolRSF

## Coast : HunterCr

## Coast : Winchuck

## Coast : ChetcoMain

## Coast : PistolRSF

## HunterCr : Winchuck

## HunterCr : ChetcoMain

## HunterCr : PistolRSF

## Winchuck : ChetcoMain

## Winchuck : PistolRSF

## ChetcoMain : PistolRSF

## user system elapsed
## 42.071 0.144 42.224
Running the AMOVA by year
yearPairs <- pairwise_amova(setpop(ramdat, ~Year), ~Year/ZONE2, quiet = TRUE)
yearPairs.cc <- pairwise_amova(setpop(ramdat, ~Year), ~Year/ZONE2,
quiet = TRUE, clonecorrect = TRUE)
set.seed(9001)
system.time(yearPairs.t <- pairwise_amova_test(yearPairs, nrepet = 999))
## 2001 : 2002

## 2001 : 2003

## 2001 : 2004

## 2001 : 2011

## 2001 : 2012

## 2001 : 2013

## 2001 : 2014

## 2002 : 2003

## 2002 : 2004

## 2002 : 2011

## 2002 : 2012

## 2002 : 2013

## 2002 : 2014

## 2003 : 2004

## 2003 : 2011

## 2003 : 2012

## 2003 : 2013

## 2003 : 2014

## 2004 : 2011

## 2004 : 2012

## 2004 : 2013

## 2004 : 2014

## 2011 : 2012

## 2011 : 2013

## 2011 : 2014

## 2012 : 2013

## 2012 : 2014

## 2013 : 2014

## user system elapsed
## 276.694 1.860 279.005
set.seed(9001)
system.time(yearPairs.cc.t <- pairwise_amova_test(yearPairs.cc, nrepet = 999))
## 2001 : 2002

## 2001 : 2003

## 2001 : 2004

## 2001 : 2011

## 2001 : 2012

## 2001 : 2013

## 2001 : 2014

## 2002 : 2003

## 2002 : 2004

## 2002 : 2011

## 2002 : 2012

## 2002 : 2013

## 2002 : 2014

## 2003 : 2004

## 2003 : 2011

## 2003 : 2012

## 2003 : 2013

## 2003 : 2014

## 2004 : 2011

## 2004 : 2012

## 2004 : 2013

## 2004 : 2014

## 2011 : 2012

## 2011 : 2013

## 2011 : 2014

## 2012 : 2013

## 2012 : 2014

## 2013 : 2014

## user system elapsed
## 44.011 0.100 44.107
Pairwise Heatmaps
By region
names(zonePairs.t) <- names(zonePairs)
names(zonePairs.cc.t) <- names(zonePairs.cc)
parray <- array(dim = c(length(zonePairs.t[[1]]$pvalue),
length(zonePairs.t),
2),
dimnames = list(Test = zonePairs.t[[1]]$names,
Pair = names(zonePairs.t),
data = c("full", "clone-censored"))
)
parray[, , "full"] <- vapply(zonePairs.t, getpval,
zonePairs.t[[1]]$pvalue)
parray[, , "clone-censored"] <- vapply(zonePairs.cc.t, getpval,
zonePairs.cc.t[[1]]$pvalue)
phiarray <- array(dim = c(length(getphi(zonePairs[[1]])),
length(zonePairs),
2),
dimnames = list(Test = rownames(zonePairs[[1]]$statphi),
Pair = names(zonePairs),
data = c("full", "clone-censored"))
)
phiarray[, , "full"] <- vapply(zonePairs, getphi,
getphi(zonePairs[[1]]))
phiarray[, , "clone-censored"] <- vapply(zonePairs.cc, getphi,
getphi(zonePairs.cc[[1]]))
p.df <- tbl_df(melt(parray[3:4, , ])) %>%
separate(Pair, c("PairOne", "PairTwo"), sep = " : ")
# Fixing the data so it appaears as an upper-triangle
setpop(ramdat) <- ~ZONE2
p.df$PairOne <- refactor(p.df$PairOne, pop(ramdat))
p.df$PairTwo <- refactor(p.df$PairTwo, pop(ramdat))
# p.df$PairOne <- unique(pop(ramdat))[match(p.df$PairOne, unique(pop(ramdat)))]
# p.df$PairTwo <- unique(pop(ramdat))[match(p.df$PairTwo, unique(pop(ramdat)))]
levels(p.df$Test) <- paste("variations between", c("year within region", "regions"))
# Plotting
ggplot(p.df, aes(x = PairTwo, y = PairOne, fill = -log(value))) +
geom_tile() +
geom_text(aes(label = value, color = value)) +
facet_grid(data ~ Test) +
scale_fill_continuous(low = "black", high = "orange", name = "p value\n",
breaks = c(0.693, 3, 4.6, 6.7),
labels = round(exp(-c(0.693, 3, 4.6, 6.8)), 3)
) +
scale_color_continuous(low = "black", high = "orange", guide = "none") +
theme_classic() + myTheme + labs(list(x = NULL, y = NULL))
## Warning: Removed 1 rows containing missing values (geom_text).
## Warning: Removed 1 rows containing missing values (geom_text).
## Warning: Removed 1 rows containing missing values (geom_text).
## Warning: Removed 1 rows containing missing values (geom_text).

# Creating a matrix to store Phi values
phi <- matrix(nrow = nlevels(p.df$PairOne), ncol = nlevels(p.df$PairOne),
dimnames = list(levels(p.df$PairOne), levels(p.df$PairOne)))
p <- phi
phi[lower.tri(phi)] <- phiarray[4, , 2] # Clone corrected (upper triangle)
phi <- t(phi) # Transposition to final form
phi[lower.tri(phi)] <- phiarray[4, , 1] # Raw data (lower triangle)
phi
## JHallCr NFChetHigh Coast HunterCr Winchuck ChetcoMain
## JHallCr NA 0.01399 0.041699 0.08363 0.03166 0.017933
## NFChetHigh 0.02571 NA 0.024833 0.07175 0.02982 0.018233
## Coast 0.06324 0.04996 NA 0.04240 0.05376 0.002047
## HunterCr 0.18133 0.14962 0.082444 NA 0.11903 0.019447
## Winchuck 0.06174 0.06309 0.071333 0.18991 NA 0.026352
## ChetcoMain 0.08349 0.08168 0.002977 -0.02980 0.07668 NA
## PistolRSF -0.01156 0.01162 0.035642 NaN 0.03453 -0.061292
## PistolRSF
## JHallCr -0.031101
## NFChetHigh -0.024162
## Coast 0.014868
## HunterCr NaN
## Winchuck -0.003834
## ChetcoMain -0.037686
## PistolRSF NA
p[lower.tri(p)] <- parray[4, , 2]
p <- t(p)
p[lower.tri(p)] <- parray[4, , 1]
p
## JHallCr NFChetHigh Coast HunterCr Winchuck ChetcoMain PistolRSF
## JHallCr NA 0.963 1.000 1.000 1.000 0.977 0.001
## NFChetHigh 0.996 NA 0.986 1.000 0.934 1.000 0.205
## Coast 1.000 0.977 NA 0.866 1.000 0.680 0.730
## HunterCr 1.000 1.000 1.000 NA 1.000 1.000 NA
## Winchuck 1.000 0.938 1.000 1.000 NA 1.000 0.667
## ChetcoMain 1.000 0.866 0.710 0.334 1.000 NA 0.341
## PistolRSF 0.001 0.809 0.720 NA 1.000 0.353 NA
write.table(phi, file = "phi_by_zone.csv", sep = ",", row.names = TRUE,
col.names = NA)
By year
names(yearPairs.t) <- names(yearPairs)
names(yearPairs.cc.t) <- names(yearPairs.cc)
parray_year <- array(dim = c(length(yearPairs.t[[1]]$pvalue),
length(yearPairs.t),
2),
dimnames = list(Test = yearPairs.t[[1]]$names,
Pair = names(yearPairs.t),
data = c("full", "clone-censored"))
)
parray_year[, , "full"] <- vapply(yearPairs.t, getpval,
yearPairs.t[[1]]$pvalue)
parray_year[, , "clone-censored"] <- vapply(yearPairs.cc.t, getpval,
yearPairs.cc.t[[1]]$pvalue)
phiarray_year <- array(dim = c(length(getphi(yearPairs[[1]])),
length(yearPairs),
2),
dimnames = list(Test = rownames(yearPairs[[1]]$statphi),
Pair = names(yearPairs),
data = c("full", "clone-censored"))
)
phiarray_year[, , "full"] <- vapply(yearPairs, getphi,
getphi(yearPairs[[1]]))
phiarray_year[, , "clone-censored"] <- vapply(yearPairs.cc, getphi,
getphi(yearPairs.cc[[1]]))
py.df <- tbl_df(melt(parray_year[3:4, , ])) %>%
separate(Pair, c("PairOne", "PairTwo"), sep = " : ")
# Fixing the data so it appaears as an upper-triangle
setpop(ramdat) <- ~Year
py.df$PairOne <- refactor(py.df$PairOne, pop(ramdat))
py.df$PairTwo <- refactor(py.df$PairTwo, pop(ramdat))
# py.df$PairOne <- unique(pop(ramdat))[match(py.df$PairOne, unique(pop(ramdat)))]
# py.df$PairTwo <- unique(pop(ramdat))[match(py.df$PairTwo, unique(pop(ramdat)))]
levels(py.df$Test) <- paste("variations between", c("regions within years", "years"))
# Plotting
ggplot(py.df, aes(x = PairTwo, y = PairOne, fill = -log(value))) +
geom_tile() +
geom_text(aes(label = value, color = value)) +
facet_grid(data ~ Test) +
scale_fill_continuous(low = "black", high = "orange", name = "p value\n",
breaks = c(0.693, 3, 4.6, 6.7),
labels = round(exp(-c(0.693, 3, 4.6, 6.8)), 3)
) +
scale_color_continuous(low = "black", high = "orange", guide = "none") +
theme_classic() + myTheme + labs(list(x = NULL, y = NULL))
## Warning: Removed 30 rows containing missing values (geom_text).
## Warning: Removed 30 rows containing missing values (geom_text).
## Warning: Removed 30 rows containing missing values (geom_text).
## Warning: Removed 30 rows containing missing values (geom_text).

phi_year <- matrix(nrow = nlevels(py.df$PairOne), ncol = nlevels(py.df$PairOne),
dimnames = list(levels(py.df$PairOne), levels(py.df$PairOne)))
phi_year[lower.tri(phi_year)] <- phiarray_year[4, , 2]
phi_year <- t(phi_year)
phi_year[lower.tri(phi_year)] <- phiarray_year[4, , 1]
phi_year
## 2001 2002 2003 2004 2005 2006 2010 2011
## 2001 NA NaN 0.041386 NaN 0 0 0 -0.02217
## 2002 NaN NA 0.033406 NaN 0 0 0 -0.01252
## 2003 -0.018524 -0.01675 NA 0.043358 0 0 0 0.02827
## 2004 NaN NaN -0.004773 NA 0 0 0 0.01991
## 2005 0.000000 0.00000 0.000000 0.000000 NA 0 0 0.00000
## 2006 0.000000 0.00000 0.000000 0.000000 0 NA 0 0.00000
## 2010 0.000000 0.00000 0.000000 0.000000 0 0 NA 0.00000
## 2011 0.031110 0.02759 0.075251 0.052300 0 0 0 NA
## 2012 -0.016712 -0.01932 -0.002365 0.018806 0 0 0 0.08346
## 2013 -0.013666 -0.02135 -0.015477 0.003758 0 0 0 0.09174
## 2014 -0.005429 -0.01381 -0.004931 0.011101 0 0 0 0.06162
## 2012 2013 2014
## 2001 -0.0182409 1.669e-03 0.0059683
## 2002 -0.0003216 -2.695e-05 -0.0042011
## 2003 0.0175266 1.206e-02 0.0146585
## 2004 0.0405752 4.187e-02 0.0311115
## 2005 0.0000000 0.000e+00 0.0000000
## 2006 0.0000000 0.000e+00 0.0000000
## 2010 0.0000000 0.000e+00 0.0000000
## 2011 0.0161974 1.903e-02 -0.0006563
## 2012 NA -1.181e-02 -0.0059601
## 2013 -0.0041298 NA -0.0034361
## 2014 0.0074928 -1.831e-02 NA
write.table(phi_year, file = "phi_by_year.csv", sep = ",", row.names = TRUE,
col.names = NA)